function [ n, result ] = Adam_Moulton( x, h, y0, funct)
%ADAM_MOULTON Summary of this function goes here
%   Detailed explanation goes here

tempx = x(1);
result = 0;

yr0 = y0;

if (nargin(funct) == 2)
    tempx_aksen = tempx + h;
    [n yr1] = Runge_Kutta([tempx tempx_aksen],h,y0,funct,4);
    tempx_aksen = tempx_aksen + h;
    [n yr2] = Runge_Kutta([tempx tempx_aksen],h,y0,funct,4);
    tempx_aksen = tempx_aksen + h;
    [n yr3] = Runge_Kutta([tempx tempx_aksen],h,y0,funct,4);
    n = ((x(2) - x(1)) / h);
    tempx = tempx + 3 * h;
    for r = 4:n
        yr = yr3 + h*(-9*funct(tempx-3*h,yr0)+37*funct(tempx-2*h,yr1)-59*funct(tempx-h,yr2)+55*funct(tempx,yr3))/24;
        tempy = yr;
        yr = yr3 + h*(funct(tempx-2*h,yr1)-5*funct(tempx-h,yr2)+19*funct(tempx,yr3)+9*funct(tempx+h,yr))/24;
        %display([tempx tempy yr]);
        yr0 = yr1; yr1 = yr2; yr2 = yr3; yr3 = yr;
        tempx = tempx + h;
    end
    result = yr;
% Fungsi 3 variabel
% ------------------------------------------------------------%
elseif (nargin(funct) == 3)
    yr0 = y0(1);
    zr0 = y0(2);
    
    tempx_aksen = tempx + h;
    [n tempy] = Runge_Kutta([tempx tempx_aksen],h,y0,funct,4);
    [n tempz] = Runge_Kutta([tempx tempx_aksen],h,y0,funct,4);
    yr1 = tempy(1);
    zr1 = tempz(2);
    tempx_aksen = tempx_aksen + h;
    [n tempy] = Runge_Kutta([tempx tempx_aksen],h,y0,funct,4);
    [n tempz] = Runge_Kutta([tempx tempx_aksen],h,y0,funct,4);
    yr2 = tempy(1);
    zr2 = tempz(2);
    tempx_aksen = tempx_aksen + h;
    [n tempy] = Runge_Kutta([tempx tempx_aksen],h,y0,funct,4);
    [n tempz] = Runge_Kutta([tempx tempx_aksen],h,y0,funct,4);
    yr3 = tempy(1);
    zr3 = tempz(2);
    n = ((x(2) - x(1)) / h);
    tempx = tempx + 3 * h;
    for r = 4:n
        yr = yr3 + h*(-9*funct(tempx-3*h,yr0,zr0)+37*funct(tempx-2*h,yr1,zr1)-59*funct(tempx-h,yr2,zr2)+55*funct(tempx,yr3,zr3))/24;
        zr = zr3 + h*(-9*funct(tempx-3*h,yr0,zr0)+37*funct(tempx-2*h,yr1,zr1)-59*funct(tempx-h,yr2,zr2)+55*funct(tempx,yr3,zr3))/24;
        %tempy = yr;
        %tempz = zr;
        yr = yr3 + h*(funct(tempx-2*h,yr1,zr1)-5*funct(tempx-h,yr2,zr2)+19*funct(tempx,yr3,zr3)+9*funct(tempx+h,yr,zr))/24;
        zr = zr3 + h*(funct(tempx-2*h,yr1,zr1)-5*funct(tempx-h,yr2,zr2)+19*funct(tempx,yr3,zr3)+9*funct(tempx+h,yr,zr))/24;
        %display([tempx tempy yr]);
        yr0 = yr1; yr1 = yr2; yr2 = yr3; yr3 = yr;
        zr0 = zr1; zr1 = zr2; zr2 = zr3; zr3 = zr;
        tempx = tempx + h;
    end
    result(1) = yr;
    result(2) = zr;    
end

end
